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ABSTRACT 

^ ' We present a method for generating virtual observations from smoothed-particle- 

hydrodynamics (SPH) simulations. This method includes stellar population synthesis 
models and the reprocessing of starlight by dust to produce realistic galaxy images. 
We apply this method and simulate the merging of two identical giant Sa galaxies 
(Mdisc = 1O U M , M S pheroid = 2.5 x 10 10 M©). The merger remnant is an elliptical 
; galaxy (M sp h e roid — 1-3 x 10 n M Q , Mdi sc — 7.4 x 1O 1O M ). The merger concentrates 

. the gas content of the two galaxies into the nuclear region. The gas that flows into the 

nuclear region refuels the central black holes of the merging galaxies. We follow the 
refuelling of the black holes during the merger semi-analytically. 
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In the simulation presented in this article, the black holes grow from 3 x 10 7 M 



to 1.8 x 10 8 M o , with a peak AGN luminosity of M B ~ -23.7. We study how the 
\ morphological and spectral properties of the system evolve during the merger and 

work out the predictions of this scenario for the properties of host galaxies during the 
active phase. The peak of AGN activity coincides with the merging of the two galactic 
nuclei and occurs at a stage when the remnant looks like a lenticular galaxy. The 
simulation predicts the formation of a circumnuclear starburst ring/dusty torus with 
/\ \ an opening angle of 30° — 40° and made of clouds with nn = 10 24 cm -2 . The average 

optical depth of the torus is quite high, but the obscuring medium is patchy, so that 
there still exist lines of sight where the AGN is visible in a nearly edge-on view. For 
the same reason, there are lines of sight where the AGN is completely obscured in the 
face-on view. 

Key words: galaxies: formation, interactions, active, nuclei - quasars: general - 
methods: N-body simulations 



1 INTRODUCTION 

iToomre Sz Toomrel (|l972f) suggested that mergers can con- 
vert spiral galaxies into elliptical ones. They also remarked 
that galaxy interactions can determine strong bar instabili- 
ties, trigger starburst s and feed gas to active galactic nuclei 
(AGN). This paradigm, which links the formation of ellip- 
tical galaxies and AGN to mergers, has derived its strength 
from its success in smoothed-particle- h ydrodynamics (SPH) 
simulations of galaxy merger s fi.e. |Bames^^Hernauist 
I l99d: iMihos fc HernauistJI 19941 : iBarnes fc HernauistJ llQQlt 
Mihos & Hernauist 1996) and in semi- analytic mod- 
els of galaxy formation fi.e. iKauffmann et alJ [l993: 
Kauffman nfc Charlotl 1998 ; Kauffmann & Haehnelt 2000; 
ICattaned 1200 J) . The discovery of ultra- luminous infrared 



galaxies (ULIRGs) with the infra-red astronomical satel- 
lite IRAS provided further support for the merger scenario. 
These galaxies release most of their power in the IR because 
most of their light is reprocessed by dust. A significant frac- 
tion of these objects host AGN and many ar e observed to be 
part of interacting or merging systems (i.e. iBushouse et alJ 

Nevertheless, from an observational point of view, the 
link of AGN with galaxy interactions is highly controver- 
sial. Since the earliest detections of quasar host galax- 
ies with the Canada-France-Hawaii Telescope, there has 
been a claim that > 30% quasar hosts are currently in - 
teracting with a companion fautchings fc CampbeJll98ab . 
IXilouris Sz PapadakiJ (|2002f) fitted de Vaucouleurs profiles 
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to HST images of faint AGN and normal galaxies from 
the local Universe. After subtracting their fits from the 
data, they concluded that all active galaxies show significant 
structure in their inner 100 pc and 1 kpc regions, contrary 
to quiescent early- type galaxies, which show no structure at 
all. However, other groups find that the isophotal profiles of 
quasar hosts are consistent with the de Vaucouleurs law and 
that their colours are those of old, passively evolving stellar 
populations, while interactions are no more co mmon than in 
norm al galaxies with similar characteristics (|Dunlop et alJ 

feoofl . 

There is now evidence that AGN seen through their 
intense optical/UV continuum (big blue bump) and broad 
emission line spectra (Type I quasars) are just the tip of 
the iceberg of the AGN phenomenon (iBrandt et al J 1 .2001: 
iBarger et all 120021 : llJeda et alJbooA ISazonov et al.1 120041^ 
at lea st in t erms of the nu mber of objects. | Yu &: Tremaine 
(I2002T) used ISoltanl dl982T) and IChokshi fc Turner! (ll9<dV s 
method to infer the cosmic density of supermassive black 
holes due to optical quasars from the quasar luminosity func- 
tion. They inferred a value of 2.1 x 10 5 M©Mpc~ 3 , which 
they compared with their best estimate of the local den- 
sity of supermassive black holes, (2.5 ±0.4) x 10 5 M©Mpc -3 , 
concluding that there is limited need for obscured or ra- 
diatively inefficient accretion. However, Barger et al. (2005) 
have shown that extrapolating a fit to the quasar luminosity 
function outside the measured magnitude range can overes- 
timate the cosmic den sity due to opt i cal qu asars by a fac- 
tor of 1.75. Moreover, iMarconi et alJ (|2004T) disagree with 
the local density of su permassive black holes estimated by 
lYu fc Tremainel fe002) and prop ose a value almost twice as 
large. Unification by orientation (jAntonucci &: Millerl fl985') 
has proven a successful paradigm to explain the type I /ll 
distinction. In this paradigm, AGN are surrounded by dusty 
tori. From a pole-on view, we see the naked AGN and the 
object is classified as type I. From an equatorial view, the 
dusty torus obscures the central engine and all we can see 
is reverberated light/narrow line emission from clouds illu- 
minated by the AGN. In this case the object is classified as 
type II. 

Studying the host galaxies of type II AGN is easier be- 
cause there is not the proble m of subtracting th e light of the 
AGN. Nature does it for us. iKauffmann et al.1 (2003) anal- 
ysed the spectra of 22 623 narrow line AGN from the Sloan 
Digital Sky Survey. By using the discontinuity at 0.4 /xm as 
an indicator of recent star formation and the strength of 
the O III 0.5007 /xm emission line as an indicator of AGN 
power, they concluded that low luminosity AGN are hosted 
by normal early- type galaxies. The hosts of bright AGN are 
galaxies with the high surface densities of early-type galax- 
ies, but bluer colours indicating recent star formation. They 
claim that this star formation is not concentrated to the nu- 
clear region. This fi nding appears to disagree with that of 
iDunlop et alJ (|2003fl and questions the paradigm of unifica- 
tion by orientation. I.e., are the hosts of broad and narrow 
line AGN drawn from the same population? If starbursts on 
the scale of the host galaxy contribute to obscure the central 
engine and the broad line region, then type I samples may 
be biased toward AGN in late- stage or gas-poor mergers. 

Here we use numerical simulations i) to work out the 
expected properties of quasar hosts in the merging scenario 
and ii) to understand the difficulties and selection effects 



that are present when we try to derive properties of quasar 
hosts from observations. For this purpose, we run hydro- 
dynamic N-body simulations of galaxy mergers, where we 
incorporate the effects of absorption by dust and a semi- 
analytic treatment of the growth of supermassive black 
holes, interfaced with visualising tools, which convert the 
outputs of the simulations into mock data (including instru- 
mental diffraction, sky noise, etc.) 

The hydrodynamic simulations use the three- 
compone nt (stars, gas and dark matter) SPH tree-code 
GalMer, (ICombes fc MelchiodEooi described in Section 2. 
The modelling of the stellar spectral energy distribution 
(SED) and of the reprocessing of star light by dust is 
performed after the simulations, through an interface that 
manipulates the simulations' outputs (Section 3). The same 
is true for the growth of the central black hole, which is 
computed semi-analytically from the rate at which gas falls 
into the central 100 pc. 

Section 5 presents the results of our analysis of the out- 
puts of the simulation. In this analysis, we concentrate on 
three questions: i) the relation between the global star for- 
mation history of the host galaxy and star formation/black 
hole accretion in the galactic nucleus, ii) the relation be- 
tween the spectro-morphological evolution of the host galaxy 
and the growth of the central black hole, and iii) the geom- 
etry, the optical depth and the covering factors of the ab- 
sorbing materials. In Section 6, we summarise the main new 
results of this work. 

This paper is part of a project to build a library of 
simulated galaxy mergers. Here, we only consider one simu- 
lation where we merge two identical Sa galaxies. In a future 
publication, we shall discuss how the results of this article 
depend on the merging parameters and on the properties of 
the merging galaxies. 

2 THE HYDRODYNAMIC SIMULATIONS 
2.1 The GalMer SPH code 

Gravitational hydrodynamics is simulated with the GalMer 
SPH tree-code, which is a development of the code used by 
ICombes &; Melchiorl |2002). In GalMer, galaxies are com- 
posed of three kinds of particles: stellar particles, dark mat- 
ter particles and hybrid particles. At the beginning of a 
simulation, a hybrid particle is made entirely of gas. We 
assume that the gas is isothermal and cools very efficiently. 
We calculate pressure forces for an isothermal gas with a 
temperature of 10,000 K. This is approximately the temper- 
ature where there is st rong change in the cooling curves 
(Blu menthal et al.l fl984V The adopted star formation rate 
is given by the modified Schmidt law: 

p* = 8.33 x 10 7 M kpc" 3 Myr _1 x 

X (lO 9 M0 S kpc- 3 ) (lOOkms- 1 ) ' (1) 

where p g is the density of the gas and a is the local ve- 
locity dispersion of the hybrid particles. We take 77 = 1.5 
for regions where the divergence is negative and the gas is 
shocked, while we take 77 = for regions where the gas is 
expanding. When the gas content of a hybrid particle drops 
below 5%, the hybrid particle is turned into a stellar particle 
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Table 1. Initial galaxy parameters in the SPH simulation 



Component 


M (Mq) 




a (kpc) 


b (kpc) 


DM halo 


1.25 • 


10 11 


40,000 





50 


Stellar disc 


1.00 • 


10 11 


48,000 


4 


0.5 


Gas disc 


1.00 • 


10 10 


20,000 


5 


0.2 


Bulge 


2.50 • 


10 10 


12,000 





2 



and its gas content is spread among the neighbouring hybrid 
particles. It is useful to think of a hybrid particle as a model 
for a giant molecular cloud, which evolves into a star cluster 
once all the gas has been locked into stars. 

2.2 The initial conditions for the merging galaxies 

In this first paper, we concentrate on the results of one 
simulation, in which we merge two identical galaxies. Each 
of the two galaxies is made of four components: the 
dark matter halo (dark matter particles), the stellar disc 
(stellar particles), the gas disc (hybrid particles) and the 
bulge (stellar particles). Each of the four components is 
a random realisation of a Miyamoto- Nagai density profile 
(iMivamoto fc NagaJl97,fl . which, in cylindrical coordinates 
r, 0, z, takes the form: 

, , = b*M ar * + [a + 3(, 2 + b^][a + (. 2 + fr 2 ) 1 / 2 ] 2 

P{ ' ' 47T {r 2 + [a+(z 2 + 6 2 ) 1 /2]2}5/2^2 + 6 2)3/2 A 

M is the total mass of the component, while a and b are 
two characteristic length- scales. The Miyamoto- Nagai den- 
sity distribution tends to a Plummer sphere with radius b 
in the limit a — > and to a Kuzmin disc with radius a in 
the limit b — > 0. These dens ity distributions are discu ssed at 
length in the textbook bv lBinnev Trem aine ( 1987). Table 
1 gives the mass M, the number of particles iV p and the 
characteristic scale-lengths a and b for the each of the four 
components. The masses and the numbers of particles that 
are given in Table 1 are the values for one galaxy. The total 
number of particles is 240,000. 

2.3 The SPH outputs 

We save the outputs of the SPH simulations at intervals of 
25 Myr. The five snapshots in Figs. 1-2 convey an idea of the 
merging dynamics. Soon after the simulation has started, the 
discs become unstable in the bar mode. The bars in the gas 
discs are stronger than the bars in the stellar discs and are 
accompanied by more violent inflows along the bar direction. 

The galaxy labelled as galaxy 1 is the one that is above 
in the 100 Myr snapshot and below in the 200 Myr snapshot 
(Fig. 1). The two galaxies are identical, but it is to this one 
that we refer when we give information pertaining to one 
galaxy only (zoom, central star formation, etc.). 



3 FROM SPH SIMULATIONS TO MOCK 
ASTRONOMICAL DATA 

3.1 Modelling stellar colours with PEGASE 

Given an initial stellar mass function (for us, a Salpeter 
IMF), the colours of a stellar population are a function of 



how old its stars are. We can reconstruct the colours of the 
hybrid particles by tracking how the gas content of each par- 
ticle varies at each 25 Myr. In this way we can decompose 
the stars of a hybrid particles into those which formed in the 
last 25 Myr, those which have an age 25 Myr < £* < 50 Myr, 
and so on. If we assume that the star formation rate is 
uniform over intervals of 25 Myr, then we can reconstruct 
the colours of the hybrid p article with the PEGASE stella r 
spectral evolution model (|Fioc &; Rocca-Volmerange]ll997f) . 
Since we only saved snapshots at intervals of 25 Myr, we lack 
the time resolution to see nebular lines. This is a purely tech- 
nical limitation dictated by the memory of the workstation 
and is easily overcome. Moreover, it has no consequence for 
the results of this paper. 

More difficult is to determine the colours of the stellar 
populations that preexisted the beginning of the simulation. 
For bulge stars we just assume a uniform age of 8 Gyr at the 
start. For the ste llar disc we follow the method described by 
ISegalovitzl (|l975f) . In his model the gas disc was assembled 
10 Gyr ago with no further growth from later gas cooling 
onto the disc. The gas is in differential rotation and, owing 
to secular instabilities, the disc very rapidly develops spiral 
arms. Each time a gas cloud passes through a spiral arm, 
a fraction p of the gas goes into stars and is not returned 
to the interstellar medium again. In this model the gas sur- 
face density evolves as H ga s(r, t) = Y, sas (r, 0) exp(— /3t). Here 
(3 = (27r) _1 n(Q— £l P )p, where n is the number of spiral arms, 
Q(r) is the rotation angular velocity required to support the 
disc and Q p is the p attern angular velocity of the spiral arms. 
ISearle et alJ |l973) determined the colours for a stellar pop- 
ulation the formation of which started 10 Gyr ago and con- 
tinued from then onwa rd at an expo nentially decaying rate 
on a time-scale /3 _1 . ISegalovitzl (|l975l) used their determi- 
nation of B — V as a function of (3 together with the data 
on B — V as a function of galactocentric radius then avail- 
able for M81. In this way he was able to reconstruct /3(r). 
He found a good fit to the data for p ~ 0.2 and Q p C Q. 
Here we use his model with n — 2, fi p = and p = 0.2 
to compute /3(r) for the discs of our model galaxies. Q(r) is 
given by the analytic formula for a M iyamoto- Nagai density 
distribution (jMivamoto &; Nagailll975fi . This scheme assigns 
an r-dependent age to each of the particles composing the 
stellar disc. PEGASE is then used to calculate the magni- 
tudes of these particles. After 10 Gyr, we expect a higher gas 
fraction in the outer parts of the disc, where matter has had 
time to complete fewer galactic rotations. We model this ef- 
fect by attributing two different scale-radii to the gas disc 
and the stellar disc (Table 1). 

Reality is more complicated because: i) the transforma- 
tion of gas into stars is accompanied by an evolution in the 
dynamic and kinematic structure of the disc, ii) star forma- 
tion is inevitably accompanied by stellar feedback, and iii) 
the outskirts of the disc will accrete fresh gas from the sur- 
rounding environment. However, the purpose of the model 
described here is simply to generate plausible initial condi- 
tions for the colours of the merging galaxies, and this simple 
model accomplishes this task quite well. 

3.2 Dust absorption/emission 

The reprocessing of star light by dust is c omputed with the 
STARDUST model (iDevriendt et allllQQflh . The optical ex- 
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Figure 1. A sequence of five snapshots from the GalMer SPH simulation. The time elapsed from the beginning of the simulation 
printed on each snapshot. The units on the axis are kpc. Black is dark matter, green is old disc stars, red is bulge stars and yellow 
gas/young stars. 
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Figure 2. The same as in Fig. 1 from an orthogonal viewing angle. 
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tinction scales linearly with the column of neutral hydro- 
gen on the line of sight and with the metallicity of the ob- 
scuring material. For a neutral hydrogen column density on 
the line of sight of nu = 10 20 cm -2 and solar abundances, 
the R-, V- and B-band extinctions are Ar c± 0.040 mag, 
Ay — 0.049 mag and Ab — 0.062 mag. The chemical evo- 
lution of the interstellar medium is already included in the 
hydrodynamic code, while nu id determined in the phase of 
post-processing, since it depends on the viewing angle. All 
hybrid particles in front of a hybrid or stellar particle con- 
tribute towards absorbing the light of that particle. More- 
over, in the case of hybrid particles, there is also a contribu- 
tion from self-absorption. 

The main difficulty derives from the fact that the 
SPH simulation cannot resolve the gas clouds that pro- 
duce the absorption. Hybrid particles correspond to gas 
clouds of 10 4 — 5 x 1O 5 M , which is the appropriate mass 
range for large molecular clouds. If the absorption of light 
was uniform over a cloud's entire cross section, then a 
cloud radius of r c i ou d would correspond to nu ~ 1.3 x 
10 21 (M cloud /10 5 M o )(r c i oud /50pc)- 2 ( assuming that 80% of 
the cloud mass M c i ou d is composed of hydrogen). However, 
in reality, molecular clouds do not absorb light uniformly be- 
cause the cold gas has a fractal distribution. Large clouds are 
composed of smaller clouds interspersed with voids. Further- 
more, the concept itself of associating hybrid particles with 
large molecular clouds is qualitatively useful but is danger- 
ous when used quantitatively because only smoothed quan- 
tities can be relied to represent physical properties. 

The more sensible approach is to compute a smoothed 
< nn > for the gas in front of each particle and then to use 
this < 71h > together with assumptions about the physical 
properties of the obscuring clouds in order to generate a 
Monte Carlo nu . It is this second nu that we eventually use 
to calculate the optical extinction. To pass from < rin > 
to a probability distribution for nu we need to make some 
hypothesis about the column densities of the clouds that 
produce the extinction. In reality, the distribution of cloud 
properties is very broad, but here we simply assume that all 
clouds have the same column density n^. If < nu ><^ n^ 1 , 
then in most cases there will be no cloud or one cloud at most 
on the line of sight. The probability of the second occurrence 
is therefore ~< un > /n^, If < >^> n^, then typically 
there will be a number of ^< nu > /n%_ clouds on the 
line of sight, and we can expect deviations from the average 
extinction to be small if this number is large. 

STARDUST is used to compute the bolometric lumi- 
nosity absorbed over the entire optical spectrum. This power 
is returned in the infrared as thermal dust emission with 
a modified blackbody spectrum. A modified blackbody de- 
scribes the physical reality better than the Planck formula 
because the dust is transparent to photons at wavelengths 
^> 200 /im, so that radiation in that part of the spectrum is 
no longer thermalised. 

When we compute broad-band infrared spectra for the 
whole galaxy, we do not consider absorption and emission 
by individual clouds. Instead, we compute a smoothed nn 
map, where one pixel corresponds to 100 pc x 100 pc. We 
also compute ± Az* and z gas ± Az gas at each pixel for 
the distribution of stellar and hybrid particles along the line 
of sight (z is the Cartesian coordinate in the direction of 
the line of sight, not to be confused with the metallicity Z). 



We use the results of these calculations to determine the 
surface brightness that is absorbed and reemitted. This is 
how we calculate the infrared emission from dust heated by 
star light. 

We also calculate infrared spectra when the AGN con- 
tributes towards heating the dust. We shall see later on that 
in our simulation ~ 80% of the lines of sight looking towards 
the central engine are obscured. The obscured lines of sight 
have optical depth r ^> 1 (see the nu column in Table 2). 
The absorbed power is therefore an important fraction of the 
AGN's bolometric luminosity. Therefore, we can make the 
approximation that there is an optically thick central region, 
which absorbs ~ 80% of the AGN power and reemits it with 
the modified black body spectrum described above. The ra- 
dius of the self-absorbed central region gives the character- 
istic radius of the surface that reemits the absorbed AGN 
light. We compute this radius to estimate the temperature 
of the reemitted radiation. 



3.2.1 Creating mock images 

The software package SkyMaker ( |http://t erap ix.ia p.fr /soft \ 
is used to convert the projected coordinates and extinction- 
corrected magnitudes of the simulated particles into mock 
images in the Johnson B, Johnson V and Cousins R bands. 
These images are generated as files in the .fits format. There- 
fore, it is possible to perform on them all the manipulations 
that could be performed on real data. SkyMaker can model 
instrumental diffraction, photon noise and atmospheric blur- 
ring. Since it is possible to put the galaxies at an arbitrary 
distance, this technology can be used not only to explore the- 
oretical scenarios for galaxies and AGN, but also to assist 
the development of new observing proposal by determining 
in advance which features would be observable at a given 
flux limit or under certain sky conditions (in the case of 
ground-based telescopes) . 

Fig. 3 shows a true colour image of one of the two galax- 
ies before the merger takes place. This snapshot is generated 
by using the RVB colour system to superimpose images in 
the B, V and R bands. This image shows the dust lane 
cutting through the galactic disc. 



4 THE BLACK HOLE 
4.1 The approach 

We assume that the two galaxies start with central black 
holes of 3xlO 7 M . This mass is inferred from the black hole 
mass to bulge mass mass relation in lHaring Sz Rixl (2004). 

The mass resolution of the simulation is not sufficient 
to compute the orbits of the black holes correctly because, 
while in reality the black hole mass is much larger than the 
stellar mass, so that the black holes sink to the centre due 
to dynamical friction, that is not the case in our simulation. 
Therefore, we impose as a constraint that the black hole is 
in the galactic nucleus at all time and calculate the accre- 
tion rate from the fuel supply that is present in the nuclear 
region. This approach relies on a correct identification the 
galactic centre. 
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usually gives a good identification of the galactic nucleus. 
Instead, by starting from the leaves with the highest den- 
sity one often misses the galactic nucleus and finds molec- 
ular clouds in the spiral arms. One can look at Fig. 1 and 
see that the centres identified by this automatic procedure, 
marked by little black crosses, are, as they should be, at the 
centres of the two galaxies. Identifications are made secure 
by the fact that the galactic nucleus is not just a name for 
the central region, but stands out as a well defined dynam- 
ical entity (see the right hand column of Fig. 5). Moreover, 
we have verified that the particles which are in the galactic 
nucleus (identified as the region at a radius of r < 100 pc) 
at the beginning of the simulation, remain in the galactic 
nucleus during the entire simulation. A small fraction are 
stripped from the galactic nucleus in mini tidal tails, but 
they eventually fall back into it. 



Figure 3. A close up image of galaxy 1 viewed edge-on shows 
the central dust lane. 



4.2 Identifying the galactic centre 

The identification of density peaks in the simulation out- 
puts is done with t he substructu r e find er adaptaHOP, de- 
scribed in detail in lAubert et al.l ([2004F ) . Basically, this al- 
gorithm applies SPH smoothing to the particle distribution, 
then finds local density maxima and saddle points in this 
smoothed distribution and constructs a tree of structures 
and substructures by using the value of the density at sad- 
dle points to control connectivity between various subcom- 
ponents of the system. The local maxima correspond to the 
leaves of the tree. A given substructure is identified as the 
set of particles verifying />sph > Psaddie, where psph is the 
SPH density associated to each particle and p sa ddie is the 
value of the SPH density at the (highest) saddle point con- 
necting this substructure to its neighbour. The output of 
adaptaHOP depends on the choice of various parameters, 
namely the number Asph of particles used to perform SPH 
smoothing, the number Ahop of neighbours used to walk 
into the particle distribution in order to find local maxima 
and to establish connectivity through saddle points, and a 
parameter /p isson controlling the significance of substruc- 
tures compared to Poisson noi se. Our choice for thes e param- 
eters is the one advocated in lAubert et alJ {2004), namely 
Asph = 64, Ahop = 16 and /poisson = 4, the latter value 
meaning that only substructures at 4cr level compared to 
pure Poisson noise are considered as real. Besides these pa- 
rameters, there is as also an overall density threshold, pth, 
which is used to eliminate all particles with SPH density 
below pth prior to the analysis. That we set here to a very 
small value (in cosmological simulations, /?th controls the 
selection of dark matter haloes). 

In adaptaHOP, a galaxy is thus viewed as a tree, which 
separates in several branches, which separate into smaller 
branches, which separate into leaves, corresponding to small 
scale density peaks. This tree represents how smaller struc- 
tures are nested into larger structures through saddle points 
at a fixed time and has nothing to do with merger trees en- 
countered in cosmological simulations of galaxy formation. 
We search for density peaks in the leaves of the tree by start- 
ing from the leaves that contain the highest masses. This 



4.3 Fuelling the central black hole by orbital 
decay of large molecular clouds 

Our SPH simulation can follow the inflow of gas from a few 
kpc into the central 100 pc, but lacks the resolution to fol- 
low the inflow of gas from a 100 pc distance into the central 
black hole. Following IShlosman &; Noguchil (|l993f) . we imag- 
ine that the stream of matter from the galactic nucleus into 
the black hole is due to the orbital decay of giant molecular 
clouds, which are subject to the dynamical friction of stars 
and of lower mass clouds and which are eventually disrupted 
by tidal forces. At that point the clouds are stretched into a 
continuous accretion flow into the black hole. The formula 
for the accretion rate from the orbital decay of molecular 
cloud is: 



M. = 2xl0- 3 ^^i 

O 1 1 05 M~ 



/ 6™ 

OT 10 5 M Q 1O 9 M ~M { 



'0 PC" 



10 km 



/ V C \^ ( a * — 

:ms _1 VlOOkms" 1 / \25kms _1 / yr 



(3) 



(IShlosman Noguchi|[jl23) , where ip is the fraction of the 
total gas in the galactic nucleus M gas that is contained in 
clouds of mass > M c i ou d, p* is the stellar density in the 
galactic nucleus, a* is the stellar velocity dispersion, v& is 
the relative velocity at which the gas rotates with respect 
to the stars and v c is the rotational velocity of the gas. We 
extract all these quantities directly from the outputs of the 
SPH simulation and use them to compute M«. The two black 
holes merge when we can no longer separate the nuclei of the 
two galaxies. The reader should be warned that the precise 
value and evolution of M. (£) is sensitive to the details of how 
we calculate the average quantities mentioned above, but we 
have verified that none of our conclusions depends on these 
details. Moreover, our model for the black hole accretion rate 
tends to smooth the accretion rate over time. The real pat- 
tern of black hole accretion is more likely to be a sequence 
of bursts with variability on all time-scales. The most prag- 
matic approach to deal with this reality when necessary is 
to introduce an AGN duty cycle, so that the quasar is only 
active for a fraction of the simulation time. The accretion 
rate during the active phase is correspondingly increased to 
account for the presence of this duty cycle. 
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Figure 4. Face-on (left column) and edge-on (right column) snapshots of the merging process. From top to bottom, the snapshots were 
taken 100, 200, 300 and 350 Myr after the simulation started. The snapshots were generated with the SkyMaker software by superimposing 
B, V and R band images. Each image corresponds to a square of 80kpc side. 
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Figure 5. The left column corresponds to left column of Fig. 4 with the difference that here we only show the young stellar population 
formed after the simulation began. In these four snapshots, and in these four snapshots only, we have neglected the extinction of starlight 
by dust. The right column shows a 2kpc x 2kpc zoom of the nucleus of galaxy 1. These snapshots were derived by zooming into the 
images on the left column of Fig. 4. Young stars dominate the emission. 
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5 RESULTS 

The evolution of the system during the merger is shown 
by the sequence of snapshots in Figs. 4-5. For each of the 
four represented time-steps, which correspond to 100, 200, 
300 and 350 Myr from the beginning of the simulation, we 
show an 80kpc x 80kpc face-on image (Fig. 4, left column), 
an 80kpc x 80kpc edge-on image (Fig. 4, right column), 
an 80kpc x 80kpc face-on image, where we only show the 
unextinguished young stellar population (the light coming 
from hybrid particles, Fig. 5, left column), and a 2 kpcx2 kpc 
zoom of the nuclear region of galaxy 1 (Fig. 5, right column). 
See the end of Section 2 for an explanation of which one is 
galaxy 1. Enlarged images of the central region are taken 
from face-on images. Compare these images with Figs. 1-2 
and notice that Sky Maker reverses the original images. 

On a large scale, the two galaxies are merging in a dy- 
namical time. Their relative orbit reaches a point of closest 
approach and has only time to make a turn before the galax- 
ies merge owing to dynamical friction. On a small scale, 
inside each galaxy, the tidal forces trigger matter inflow. 
This process culminates with the merger of the two nuclei, 
which occurs about 50 Myr after the galaxies as a whole 
have merged, and ends with a stream of debris that keeps 
on falling into the nuclear region during the post-merger 
phase of dynamical relaxation. 

5.1 Star forming regions 

The most obvious result coming out of these images is that 
strong star formation is concentrated in a series of knots 
along gas spiral arms, which are much more pronounced 
and filamentary than the spiral arms in the stellar disc. The 
galactic nucleus stands out as the central one and the most 
prominent of these knots. As the simulation progresses, these 
knots, associated to the formation of massive star clusters, 
fall to the centre and merge with the galactic nucleus. 

The star formation rate in the host galaxy is a good 
tracer of the star formation rate in the nuclear region 
(Fig. 6). For this reason, we do not find evidence for a delay 
between star formation in the host galaxy and supply of fuel 
to the central region. 

The first peak of star formation in Fig. 6 is exaggerated 
by the artificial assumption that the discs are initially ax- 
isymmetric (Fig. 7). The assumed initial condition is highly 
unstable. The first peak of star formation would have been 
absent if we had started from a relaxed dynamical config- 
uration. The second peak corresponds to a physical phe- 
nomenon, the starburst triggered by the merger. During this 
starburst the star formation rate is about six times higher 
than it is in an isolated galaxy. 

5.2 The dusty torus 

The second important finding is the formation of a nuclear 
starburst ring or dusty torus at the galactic centre. The large 
number of particles in the central 300 pc combined with the 
use of a tree code enable us to obtain a higher resolution 
in the central region. We have estimated that, even under 
the most disfavourable assumptions, our resolution limit in 
the central region is smaller than 40 pc. Therefore, we are 
confident to resolve the hydrodynamics of the central 300 pc 
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Figure 6. The global star formation rate in the merging system 
(solid line) compared to the star formation rate in the nucleus of 
galaxy 1 (dotted line). 




Myr 

Figure 7. Ratio between star formation rate in the merger simu- 
lation and star formation rate for an isolated galaxy as a function 
of time. 



region, from the mean properties of which we calculate the 
black hole accretion rate semi-analytically (Eq. 3). The torus 
in Fig. 5 has an inner radius of ~ 40 pc and an outer radius 
~ 250 pc with most of the mass at ~ 130 pc. The mass- 
averaged rotation speed and the one- dimensional velocity 
dispersion within a distance r max = 250 pc from the centre 
of the remnant are v TO t — 140kms _1 and a — 63kms -1 , 
respectively. We have verified that v^ ot c± GM(r max )/w, 
where M(r max ) is the mass at r < r max and G is the gravita- 
tional constant, which means that the torus is supported by 
rotation. The ratio between the thickness and the diameter 
of the torus is consistent with the a/v TO t ratio. The forma- 
tion of a central starburst ring such as t he one described in 
this article has already been reported bv lHeller fc Shlosman 
(|l994T) . They simulated the dynamical evolution of globally 
unstable galactic discs and discussed how supernova explo- 
sions, radiation driven winds from massive stars and the 
presence of supermassive black holes affect the dynamical 
evolution of the starburst ring. Our simulation does not 
model these effects, which may play an important role in the 
evolution of the central region. We observe the formation of 
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Table 2. Virtual observations of the AGN at different viewing 
angle 
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a well defined torus after the nuclei of the two galaxies have 
merged, but after ~ 50 Myr the torus is no longer recognis- 
able in an obvious manner. At that point the structure of 
the galactic nucleus changes from a torus into a mini-spiral 
with asymmetric spiral arms as more gas and stars fall into 
the central region. Less obvious tori are also seen before the 
two nuclei merge. 

The nuclear region is very dusty. We have assumed that 
the torus is made of clouds with = 10 24 cm -2 (like the 
Orion nebula) and we have conducted virtual observations 
of the AGN 350 Myr after the simulation began. We have ob- 
served the AGN at inclination angles from pole-on (0 = 0°) 
to edge-on (0 — 90°). The azimuthal angle was selected 
randomly at each observation. In Table 2 we give the results 
of this Monte Carlo experiment showing the smoothed col- 
umn density < nu > and the extinguished AGN blue magni- 
tude Mb for each observation. Mb = —23.7 is the absolute 
blue magnitude of the unobscured AGN and Mb — 376 cor- 
responds to one cloud on the line of sight. 

The average column density of neutral hydrogen on the 
line of sight to the AGN is fairly high at all viewing angles 
and grows by an order of magnitude when we pass from a 
pole-on to an edge-on view. If the dusty torus is composed 
of clouds with high optical depth, then the AGN appears as 
a Seyfert 1 when it is observed pole-on and as a Seyfert 2 
when it is observed edge-on. Table 2 shows that there are 
a few exceptions to this rule due to presence of clouds in 
the polar regions and holes in the dusty torus. Table 2 also 
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Figure 8. Surface density profile of the 1st galaxy 100 (dashed 
line), 200 (dotted line), 350 (dashed-dotted line) and 500 (solid 
line) Myr after the simulation has started. 
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Figure 9. Surface density profile of the merger remnant 500 Myr 
after the simulation has started (crosses). The solid line fit is 
the sum of a de Vaucouleurs profile oc exp[— (r/r e ) 1//4 ] (dotted 
line) and an exponential disc oc exp(— r/ra) (dashed line). The 
fit shown on this diagram corresponds to M spheroic i = 1.34 x 
10 n M Q , M disc = 7.4 x 1O 1O M . 

shows that the metallicity of the dusty torus is about twice 
the solar value due to intense star formation activity in the 
torus. 

5.3 Quasar hosts 

Quasar host galaxies have been investigated according to 
both morphological and spectral criteria. We have analysed 
the morphology of the merger remnant by using IRAF to fit 
elliptical isophotes to R-band mock images. We compared 
the photometric profiles extracted by this procedure with de 
Vaucouleurs and exponential disc profiles. We repeated the 
same analysis by using mass surface density profiles deter- 
mined by computing the surface density in circular coronae 
around the galactic centre. Both methods gave similar re- 
sults in terms of the bulge/disc decomposition. Fig. 8 shows 
how the surface density profile evolved from t = 100 Myr to 
t = 500 Myr. We see a transition from a more exponential 
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Figure 10. Mock Hubble Space Telescope image of the system 
500 Myr after the simulation has started. The system is assumed 
to be at z = 0.27 and is observed at a viewing angle from which 
it is possible to see the central engine, hence the diffraction rings. 
The imaging model mimics the 814 jum filter on the Wide Field 
Camera. 



disc- type to a more de Vaucouleurs-type profile. In Fig. 9 
we show a spheroid plus disc decomposition of the surface 
density profile at 500 Myr. This decomposition corresponds 
to M S p h e r oid = 1-34 x 10 n M Q and M dis c = 7.4 x 1O 1O M . 
Fig. 10 shows how the merger remnant would appear in 
a Hubble Space Telescope image if it was at a redshift of 
0.27 and was photographed with the Wide Field Camera 
mounted with the 814 fim filter. The AGN is visible and has 
a luminosity of Mb — —22. We clearly see two tails of debris 
falling on the merger remnant from opposite directions. 

Fig. 11 shows how the optical and infrared spectrum 
evolves during the merger. The solid lines show the total 
spectrum while the dashed lines show the spectrum with- 
out the AGN. The optical spectrum is produced by direct 
light from stars and from the AGN. The infrared/sub- mm 
spectrum is produced by thermal dust heated by stars and 
by the AGN. The AGN contributes to the infrared and 
sub-millimetre emission independently of the viewing angle 
because most of the AGN power is absorbed and thermal- 
ized whether the c louds are on the line of sight or not (e.g. 
bazonov et al.l l2004V On the other hand, the direct optical 
emission from the AGN is totally extinguished in most cases. 

In Fig. 11, B-band corresponds to Log(A//im) ~ —0.37. 
At all times, the old stellar populations the most impor- 
tant contribution to the optical spectrum comes from the 
old stellar population preexisting the merger, not only be- 
cause most of the stellar mass is in old stars, but also be- 
cause young stars are heavily obscured by dust. However, 
as the merger progresses, clear signals of star forming ac- 
tivity emerge from the infrared and ultraviolet parts of the 
spectrum. The presence of an AGN raises the infrared flux 
substantially by increasing the temperature of the dust in 
the central starburst ring from 43 K to 63 K. 



6 DISCUSSION AND CONCLUSION 

We have developed a tool that can simulate real galaxy im- 
ages to a very high degree of realism. Our study confirms 
the powerful effect of mergers on galactic morphologies and 
the accumulation of cold gas i n the galactic centre that 
were seen in previous stu dies fi.e. lBarnes Hernauistll99ll : 
Mihos & Hernauist 199^). The gas that concentrates in the 
galact ic centre forms a s t arbur st ring, which can fuel an 
AGN. iHeller Shlosmanl (l994) have already encountered 
this configuration in SPH simulations of globally unstable 
isolated disc galaxies. They have made simulations with and 
without supernova feedback. Without feedback, gas has a 
much stronger tendency to fragment into few large clumps, 
but the effect on the final black hole mass is not significant. 
We are aware that we lack the resolution to investigate the 
distribution of gas in the starburst ring and the black hole 
accretion rate from the SPH simulation, although our reso- 
lution length is at least three times smaller than the charac- 
teristic radius of the torus. Instead, we study these aspects 
by applying a semi-analytic treatment to the outputs of the 
simulation. 

Even under conservative assumptions for the fuelling 
rate, the accretion of gas increases the black hole mass by a 
factor of > 3. Another factor of 2 is gained by merging the 
two black holes, so that the black hole starts with a mass of 
3 x 10 7 Mq and reaches a mass of - 1.8 x 10 7 Mq by the end 
of the simulation. We have simulated a relatively gas-poor 
merger. The starburst would be much stronger if the discs 
had a higher gas fraction. So would the growth of the central 
black hole. 

This increase in black hole mass is comparable to the 
increase in bulge mass from 2.5 x 1O 1O M to 1.34 x lO n M , 
corresponding to a factor of c± 5.4. Therefore, with our nor- 
malisation of the black hole accretion rate in Eq. (3), the 
black hole mass to bulge mass ratio is comparable at the 
beginning and at the end. The observations suggest that 
the scatter in the black hole mass to bulge mass relatio n 
is small (see the Introduction and also lHaring fc R ix 2004). 
That means that our normalisation of the black hole accre- 
tion rate is not far away from the appropriate value because 
it preserves the M./M sp heroid ratio. 

The star formation rate in the nuclear region closely fol- 
lows the global star formation of the host galaxy. It has two 
peaks, which coincide with the first and the second pericen- 
tre passage. Fig. 7 compares the star formation in our merger 
simulation with the star formation rate in a simulation in 
which there is only one isolated galaxy. Manifestly, the dy- 
namical instability deriving from the initial conditions for 
the discs, which are thin and axisymmetric, exaggerate the 
pace at which spiral arms form and gas fragments into knots. 
Nevertheless, we believe that the filamentary lumpiness of 
star formation evident in Fig. 5 reflects a real astrophysical 
phenomenon. Therefore, the first peak in Fig. 6 is artificial. 
Instead, the second peak is physical and corresponds to the 
actual merger. The starburst triggered by the merger is ^ 6 
times more intense than quiescent star formation. Its du- 
ration is of the order of 10 Myr. This peak coincides with 
maximum AGN activity and occurs when the merger rem- 
nant begins to look like a lenticular galaxy both colour- wise 
and morphology- wise (Fig. 4). 

Young stars are concentrated in the nuclear starburst 
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Figure 11. Global spectral energy distribution of the interacting system 100, 200, 300 and 350 Myr after the simulation has started. At 
optical wavelengths, the thick line (above) corresponds to the case in which in the AGN is visible and the thin line (below) to the case 
in which it is obscured. At infrared wavelengths, the thick line (above) corresponds to dust heated by stars and by the AGN. The thin 
line (below) corresponds to dust heated by stars only. 



ring. Therefore, it is difficult to disentangle their light from 
the AGN point-spread-function (point sources appear to 
have spikes connected by a spider web owing to instrumen- 
tal diffraction, Fig. 10). However, this is a merger with a 
low gas fraction. It is possible that, in the case of a gas- 
rich merger, star formation is no longer confined to the 
nuclear region, but dis tributed throughout the host galaxy 
dKauffmann et aljEoO^l . 
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We observe the formation of a dusty torus, which 
obscure the AGN from most viewing angles. The torus 
has as an opening angle of 30° — 40° and its proper- 
ties are consistent with expectations of unified model (see 
Kowan-KoDinsc if the torus is made of clouds with 

71h ~ 10 24 cm -2 (like the Orion nebula). Our model predicts 
that in some case the AGN is obscured even from a face-on 
viewing angle owing to the presence of clouds in the polar 
regions, while, in other cases, we can still see the AGN from 
nearly equatorial viewing- angles because there are holes in 
the dust distribution. 
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